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Abstract 

This paper introduces generalized Sobol' indices, compares strategies 
for their estimation, and makes a systematic search for efficient estimators. 
Of particular interest are contrasts, sums of squares and indices of bilinear 
form which allow a reduced number of function evaluations compared to 
alternatives. The bilinear framework includes some efficient estimators 
from Saltelli (2002) and Mauntz (2002) as well as some new estimators 
for specific variance components and mean dimensions. This paper also 
provides a bias corrected version of the estimator of Janon et al. (2012) and 
extends the bias correction to generalized Sobol' indices. Some numerical 
comparisons are given. 



1 Introduction 

Sobol' indices are certain sums of variance components in an ANOVA decom- 
position. They are used to understan d the importa n ce of various subsets of 



variables in global sensitivity analysis. ISaltelli et al.l ([20081 ) give an extensive 



introduction to Sobol' indices and variance based methods in general for investi- 
gating computer models. Linear combinations of Sobol' indices are also used to 
measure the effective dimension of functions for quasi- Monte Carlo integration. 

This article reviews Sobol' indices for a statistical audience, relating them 
to well known ideas in experimental design, particularly crossed random effects. 
Moving from physical experiments to computer experiments brings important 
changes in both the costs and goals of the analysis. In physical experiments one 
may be interested in all components of variance, or at least all of the low order 
ones. In computer experiments interest centers instead on sums of variance 
components. While the ANOVA for computer experiments is essentially the 
same as that for physical ones, the experimental designs are different. 

Of particular interest are what are known as 'fixing methods' for estimation 
of Sobol' indices. These evaluate the function at two points. Those two points 
have identical random values for some of the input components (the ones that 
are 'fixed') but have independently sampled values for the other components. 
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Sample variances and covariances of point pairs are then used to estimate the 
Sobol' indices. 

As a basic example, let / be a deterministic function on [0, l] 5 . One kind of 
Sobol' index estimate takes a form like 

Cov(/(xi, x 2 , x 3 , Xi, x 5 ),f(xi,x 2 , x 3 , Zi, z B )) (1) 

for Xj ~ U(0, 1) independently of Zj ~ U(0, 1). As we will see below, this index 
measures the sum of variance components over all subsets of the first three input 
variables. The natural design to estimate (JlJ consists of n pairs of function 
evaluations, which share the first 3 input values, and have independe nt draws in 



the la st 2 inputs. In the language of statistical experimental design (|Box et al 



19781) this corresponds to n independently generated lxlxlx2 2 1 designs. The 



first three variables are at 1 level, while the last two are a fractional factorial. 
Each replicate uses different randomly chosen levels for the five variables. 
An example of the second kind of Sobol' index is 

-Va.r(f(xi,X2,X3,Xi,x 5 ) - f(xi,x 2 ,x 3 ,Z4,,z 5 )) (2) 
= Var(/(cci, x 2 , x 3 , x±, x 5 )) - Cov(f(xi,x 2 ,x 3 , x\, x 5 ), f{x\, x 2 , x 3 , z±, z 5 )) 

The sampling design to estimate ^ is that same as that for , but the quantity 
estimated is now the sum of all variance components that involve any of the first 
three variables. The difference between ([T]) and ([2]) is that the latter includes 
interactions between the first three and last two variables while the former 
excludes them. 

The great convenience of Sobol's measures is that they can be directly es- 
timated by integration, without explicitly estimating all of the necessary in- 
teraction effects, squaring them, integrating their squares and summing those 
integrated squared estimates. Sobol' provides a kind of tomography: integrals 
of cross-products of / reveal facts about the internal structure of /. 

The goal of this paper is to exhibit the entire space of linear combinations 
of cross-moments of function evaluations with some variables fixed and others 
independently sampled. Such a linear combination is a generalized Sobol' index, 
or GSI. Then, using this space of functions, we make a systematic search for 
estimators of interpretable quantities with desirable computational or statistical 
properties. 

This systematic approach yields some new and useful estimators. Some have 
reduced cost compared to previously known ones. Some have reduced sampling 
variance. It also encompasses some earlier work. In particul ar, an e f ficient 
strategy to estimate all two factor interaction mean squares due to ISaltellil (|20f)2h 
appears as a special case. 

Section [2] introduces some notation and reviews the ANOVA of [0, l] d and 
Sobol' indices. A compact notation is necessary to avoid cumbersome expres- 
sions with many indices. Section [3] defines the generalized Sobol' indices and 
gives an expression for their value. It also defines several special classes of GSI 
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based on interpretability, computational efficiency, or statistical considerations. 
These are contrasts, squares, sums of squares and bilinear GSIs. Section[4]shows 
that many GSIs including the Sobol' index ([1} cannot be estimated by unbiased 
sums of squares. Section [S] considers estimation of a specific variance component 
for a proper subset containing k of the variables. A direct approach requires 2 k 
function evaluations per Monte Carlo sample, while a bilinear estimate reduces 
the cost to 2 L fc / 2 J -|_2 fe -L fe / 2 J . That section also introduces a bilinear estimate for 
the superset importance measure defined in Section [5] Section [5] considers some 
GSIs for high dimensional problems. It includes a contrast GSI which estimates 
the mean square dimension of a function of d variables using only d + 1 func- 
tion evaluations per Monte Carlo trial as well as some estimators of the mean 
dimension in the truncation sense. Section [7] presents a bias correction for GSIs 
that are not contrasts. Section [8] makes some comparisons among alternative 
methods and Section |H] has conclusions. 



2 Background and notation 



The analysis of variance originates with IFisher and Mackenzie! ( 1923 ) . It parti- 
tions the variance of a quantity among all non-empty subsets of factors, defined 
on a finite Cartesian grid. 

The ANOVA was generalized bv lHoeffding] (| 19481 ) to functions in L 2 [0, l] d for 
integer d ^ 1. That generalization extends the one for factorial experimental 
designs in a natural way, and c an be applie d to L 2 functions on any tensor 



product domain. For d — oo, see lOwenl (|1998l ). 



The ANOVA of L 2 [0, l] d is also attributed to ISoboll (|l969f ). For historical 
interest, we note that Sobol' used a different approach than Hocffding. He 
represented / by an expansion in a complete orthonormal basis (tensor products 
of Haar functions) and gathered together terms corresponding to each subset of 
variables. That is, where Hoeffding has an analysis, Sobol' has a synthesis. 

We use x = (xi,X2, ■ ■ ■ , Xd) to represent a typical point in [0, l] d . The set 
of indices is T> = {1,2,..., d}. We write it C v to denote a proper subset, that 
is u C v. For u C T> we use \u\ to denote the cardinality of u, and either —u 
or u c (depending on typographical clarity) to represent the complementary set 
T> — u. The expression u + v means u U v where u and v are understood to be 
disjoint. 

For u — ■ ■ ■ ,j\ u \} ^ V the point x u e [0,1]'"' has components 



(Xj 1 , Xj 2 



). The differential dx u is Y\^ eu dxj 



The ANOVA decomposition represents f(x) via 

fix) = J2 /«(*) 

uCV 

where the functions /„ are defined recursively by 



fu(x) 



I (f(x)-J2fv(x))dx 
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= / f(x)dx- u -y2f v (x). 

J[0,l] d -^ vC u 

In statistical language, u is a set of factors and f u is the corresponding effect. We 
get f u by subtracting sub-effects of f u from / and then averaging the residual 
over Xj for j u. 

From usual conventions, fz{x) = fx = f, Q ^ f(x) dx for all x e [0, l] d . 
The effect /„ only depends on Xj for those j <E u. For / e L 2 [0,l] d , these 
functions satisfy J Q f u (x)dxj — 0, when j 6 u (proved by induction on 
from which it follows that J f u (x)f v (x) dx = 0, for The ANOVA identity 

is cr 2 = J2u a u where cr 2 = J (/(a;) — /i) 2 dcc, cr 2 , = and cr 2 = J /„(a;) 2 da; is 
the variance component for m^0. 

We will need the following quantities. 

Definition 1. For integer d ^ 1, let u and v be subsets of T> — {1, . . . , d}. Then 
the sets XOR(w, v) and NXOR(u, v) are 

xor(m, u)=uU!i-nni) 
NXOR(w, v) = {u n «) (J (u c n w c ). 

These are the exclusive-or of u and i> and its complement, respectively. The 
NXOR operation satisfies the following easily verifiable properties: NXOR(u, v) = 
NXOR(w, u), and if w = NXOR(u, v) then u = NXOR(w, w). 



2.1 Sobol' indices 

This section introduces the Sobol' indices that we generalize, and mentions some 
of the methods for their estimation. 

There are various ways that one might define the importance of a variable 
Xj. The importance of variable j <G {1, . . . ,d} is due in part to c 2 ^}, but also 
due to cr 2 for other sets u with j S u. More generally, we may be interested in 
th e impor t ance of x u for a subset u of the variables. 

ISobol'l (Jl993) introduced two measures of variable subset importance, which 
we denote 

i£ = $^ ct «' and t 2 u = a l- 

We call these the lower and upper Sobol' index for the set u, respectively. The 
lower index is a total of variance components for u and all of its subsets. One 
interpretation is t 2 = Va,T(/j, u (x)) where (J, u {x) — E(/(cc) | x u ). 

The upper index counts every ANOVA component that touches the set u in 
any way. If r 2 is large then the sub set u is clearly impo rtant. If r 2 is small, then 
the subset u is not important and ISobol et al. (2007) investigate the effects of 



fixing such x u at some specific values. The second measure includes interactions 
between x v and x^ v while the first measure does not. 
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These Sobol' indices satisfy r„ ^ f u and t_u + T— u = a ■ Sobol' usually 
normalizes these quantities, yielding global sensitivity indices r^/cr 2 and r 2 /a 2 . 
In this paper we work mostly with unnor malized versio ns. 

Sobol ' s origi n al work was published inlSobol'l (Il990h before being translated 



Soboll (|1993I ). Ilshigami and Hommal (|1990l ) independently considered com- 



putation of r 2 



To estimate Sobol' indices, one pairs the point x with a hybrid point y that 
shares some but not all of the components of x. We denote the hybrid point by 
y = x u :z- u where yj = Xj for j g u and y u — Zj for j £ u. 

From the ANOVA properties one can show directly that 

/ /(x)/(x„:z_„) dec dz_„ = V" / f v (x)f(x u :z- u )dxdz- u 



[04] 2 

2 i 2 
A« +T.U- 



This is also a special case of Theorem[2]below. As a result r 2 = Cov(/(x), f(x u :z- u )) 
and so one can estimate this Sobol' index by 



1 - 

zl = -]Tf(xi)f(xi, v :Zi,- 

71 * * 



A 2 , 



(3) 



for Xi,Zi ~ U(0, l) d , where fi — (l/n)X^Li f( x i)- ^ is e ven better to use [i = 
(1/n) + f(xi tU :Zi - u ))/2 instead, as shown by I Janon et al. ( 2012 ). 

Similarly one can show that 



[0,l] d +l l 



(/(*) - f{x-u-z u )) 2 dx dz u , 



and one gets the estimate 

2 1 " 



for ! z i) ; Si i ^U(0,l) d . 

The estimate T u is unbiased for f 2 u but r 2 above is not unbiased for t 2 . It has 
a bias equal to — Var(/t). If J" |/(a;)| 4 da; < oo then this bias is asymptotically 
ne gligible, but in cases where t 2 is small , the b ias may be important. 



Mauntz ( 20021 ) and lKucherenko et al. ( 2011 ') use an estimator for t 2 derived 
as a sample version of the identity 



f{x){f{x u -z- u ) - f{z))dx dz. 



(4) 



Saltellil ([2002) also mentions this estimator. Here and below, integrals are by 
default over x E [0, l] d and/or z € [0, l] d even though some components Xj or 



Zj may not be required. An estimator based on Q with Xi,Zi ~ U(0, l) d is 
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unbiased for r„, but it requires 3n evaluations of / instead of the 2n required 
by ((3j) with either formula for fi. Proposition [3] in Section [7] gives an unbiased 
variant on the estimator of Q using only 2 function evaluations per (xi,Zi) 



pair. 



There are 2 d — 1 variance components a 2 as well as 2 d — 1 Sobol' indices 
and r„ of each type. We can recover any desired a 2 as a linear combination of 



For example a 2 ^ 2 3 | 



More generally, we have the Moebius-type relation 

^ = Ec-i) 1 "- 1 ^ . 



{1,2,3} -{1,2} -{1,3} -{2,3} ^-{1} ^-{2} ^-{3}" 



(5) 



Because / is defined on a unit cube and can be computed at any desired 
point, methods other than s imple Monte Carlo can be applied. Quasi-Monte 
Carlo (QMC) sam pling (s e e iNiederreiterl (|l992l) ) can be used instead of plain 



Monte Carlo, and ISobol'l (|200ll ) reports that QMC is more effective. For 
functions ,f that are very ex pensive, a Bayesian numerical analysis approach 
(jOaklev and O'Haganl |200J) based on a Gaussian process model for / is an 



attractive way to compute Sobol' indices. 



2.2 Related indices 

Another measure of the importance of x u is the superset importance measure 



E^ 2 



(6) 



used by iHooker ( 2004 ) to quantify the effect of dropping all interactions con- 
taining the set u of variables from a black box function. 

Sums of ANOVA components are also used in quasi-Monte Carlo sampling. 
QMC is, in general, more effective on integrands / that are dominated by their 
low order ANOVA components. Two versions of / that are equivalent in Monte 
Carlo sampling may behave quite differently in QMC. For example the basis 
used to s ample Brownian paths has been seen to affect the accuracy of QMC 
integrals (|Caflisch et all Il997t lAcworth etaD . Il997i : llmai and Tail l2002h. 



The function / has effective dimension s in the superposition sense (jCaflisch et al 



19971) . if J2\u\<8 a l > (1 - e )< j2 - Typically e = 0.01 is us ed as a default . Sim i 



larly, / has effective dimension s in the truncation sense (jCaflisch et all Il997l ) 



if E„C{1,2,..., S } °l = !{1,2 s} > (1 - £ ) ff2 - 

It is much easier to estimate the mean dimension (superposition sense) de- 
fined as ~^2 U \u\o\l a 2 than the effective dimension, because the mean dimension 
is a linear combination of variance components. The mean dimension also offers 
better resolution than effective dimension. For instance, two functions having 
identical effective dimension 2 might have different mean dimensions, say 1.03 
and 1.05. Similarly one can estimate a mean square dimension \u\ 2 / a 2 : 
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Theorem 1. 



£ii } = £« (7) 



j — l k^j u u 



Proof. This follows from Theorem 2 of Liu and Owen ( 20061) . □ 



3 Generalized SoboP indices 

Here we consider a general family of quadratic indices similar to those of Sobol'. 
The general form of these indices is 

^ ^ ^ ^ ^uv I I f (p^u -Z — u) fip^v -Z — v 

)dxdz (9) 

uCVvCV J J 

for coefficients Sl uv . If we think of f(x) as being the standard evaluation, then 
the generalized Sobol' indices © are linear combinations of all possible second 
order moments of / based on fixing two subsets, u and v, of input variables. 
A matrix representation of © will be useful below. First we introduce the 

Sobol' matrix 9 € M 2 x2 with entries 9 ut , = ff f(x u :z- u )f(x v :z- v )dxdz 
indexed by subsets u and v. Then @ is the matrix inner product tr(J7 T 0) for 
the matrix ft. Here and below, we use matrices and vectors indexed by the 2 d 
subsets of T>. The order in which subsets appear is not specified; any consistent 
ordering is acceptable. 

The Sobol' matrix is symmetric. It also satisfies Q uv = & u c v c - Theorem [2] 
gives the general form of a Sobol' matrix entry. 

Theorem 2. Let f E L 2 [0, l] d for d ^ 1, with mean /j = j f(x) dx and variance 
components for mCP. Let u,v C T>. Then the uv entry of the Sobol' matrix 
is 

= A* + 21nxor(u,u)- 

Proof. First, 

6u„ = JJ f(x u :z- u )f(x v :z- v )dxdz 

= £ E fw( x u-z-u)fw>(x v :z- v )dxdz 

wCVw'CV' J 



wCV 



E fw{x u :z- u )f w (x v :z_ v )dxdz, 



because if w ^ w' , then there is an index j G XOR(w, w'), for which either the 
integral over Xj or the integral over Zj of f w f W ' above vanishes. 
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Next, suppose that w is not a subset of NXOR(w, v). Then there is an index 
j € XOR(u, v) fl w. If j € w n it n v c , then the integral over Xj vanishes, while 
if j G w fl u (1 u c , then the integral over Zj vanishes. Therefore 

®uv ^ II fw(x u .Z— u ^f w (z v .X — v ^dxdz 

w<ZNXOR(u,v) 



E 

«iCNXOR(u,») 

2 , 2 
/i -h lNXOR(ti,n)' 



y f w (x) 2 dx 



□ 



3.1 Special GSIs 

Equation ([9]) describes a 2 2d dimensional family of linear combinations of pair- 
wise function products. There are only 2 d — 1 ANOVA components to estimate. 
Accordingly we are interested in special cases of (|9|) with desirable properties. 

A GSI is a contrast if Ylv ^ uv = Contrasts are unaffected by the 
value of the mean /i, and so they lead to unbiased estimators of linear com- 
binations of variance components. GSIs that are not contrasts contain a term 
\j? J2 U J2v ^uv an d require us to subtract an estimate fl 2 Y^v ^™ m or der 
to estimate Y^ u J2v ^™^NXOR(u,t>) as in Section [7] 

A generalized Sobol' index is a square if it takes the form 



X u f(x u :Z- u ) dxdz 



Squares and sums of squares have the advantage that they are non-negative and 
hence avoid the problems associated with negative sample variance components. 
A square GSI can be written compactly as tr(AA T 9) = A T 9A where A is a vector 
of 2 d coefficients. If X u is sparse (mostly zeros) then a square index is inexpensive 
to compute. 

A generalized Sobol' index is bilinear if it takes the form 

^2 Kf(x u :z^ u ) J ( ^2~f v f(x v :z^ v ) j dxdz. 

U V ' 

Bilinear estimates have the advantage of being rapidly computable. If there 
are ||A||o nonzero elements in A and H7H0 nonzero elements in 7 then the in- 
tegrand in a bilinear generalized Sobol' index can be computed with at most 
IItIIo + ll^llo function calls and sometimes fewer (see Section EOf even though 
it combines values from U7H0 X ||A||o function pairs. We can write the bilinear 
GSI as tr(A7 T 0) — 7 T 9A. The sum of a small number of bilinear GSIs is a low 
rank GSI. 

A GSI is simple if it is written as a linear combination of entries in just one 
row or just one column of 9, such as 



^2 X u f(x u :Z- u )f{z) dx dz. 
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It is convenient if the chosen row or column corresponds to u or v equal to 
or T>. Any linear combination ^2 u <5 M (/i 2 + r 2 ) of variance components and jj, 2 
can be written as a simple GSI taking X u = 5- u . There are computational 
advantages to some non-simple representations. 



3.2 Sample GSIs 

To estimate a GSI we take pairs (xj,z,) ~ U(0, l) 2d for i = l,...,n and 
compute tr(f2 T 0) where 



1 n 

n * — ' 



n 

i=l 



We can derive a matrix expression for the estimator by introducing the vectors 

Fi = F( Xi ,Zi) = (f{ Xi , u : Zi ,- u )) u ^ v G K 2dxl , 
for i = 1 , . . . , n and the matrix 

F = (Fi F 2 ••• f„) T el" x2i . 

The vectors Fi have covariance 8 — ^l 2 . Then 

1 ™ 1 

e = -Yf 1 fJ = — f t f, 

i=l 



and the sample GSI is 

tr(f2 T 6) = tr(9 T 0) = -Ur(F T Fft) = ^tr(F T fiF). 

3.3 Cost per (x, z) pair 

We suppose that the cost of computing a sample GSI is dominated by the 
number of function evaluations required. If the GSI requires C(il) (defined 
below) distinct function evaluations per pair (x i} Zi) for i = l,...,n, then the 
cost of the sample GSI is proportional to nC(fi). 

If the row il uv for given u and all values of v is not entirely zero then we 
need the value f(x u :z_ u ). Let 



C„.(fi) 



1, 3v<ZV with Q uv ^ 
0, else 



indicate whether f(x u :z- u ) is needed as the 'left side' of a product f(x u :z- u )f(x v :z- v ). 
The number of function evaluations required for the GSI tr(f2 T 9) is: 

c(0)= Y<( c u.(n) + c u .(n T )-c u .(n)c u .(n T )y (10) 

uCV 

We count the number of rows of Q for which f(x u :z- u ) is needed, add the 
number of columns and then subtract the number of double counted sets u. 
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4 Squares and sums of squares 



A square or sum of squares yields a nonnegative estimate. An unbiased and 
nonnegative estimate is especially valuable. When the true GSI is zero, an 
unbias ed nonnegative estimate will always return exactly zero as Fruth et al.l 
( 20121 ) remark. The Sobol' index f\ is of square form, but is not. Theorem[TJ 
lea ds to a sum of s quare s for J2 U \ u \ a u- 



Liu and Owen! (|2006f ) express the superset importance measure as a square: 



2l"l 



dx dz. 



(11) 



Fruth et al. ( 20121 ) find that a sample version of (TTT|) is the best among four 
estimators of T 2 , . 



In classical crossed mixed effects models ([Montgomervl 119981 ) every ANOVA 
expected mean square has a contribution from the highest order variance com- 
ponent, typically containing measurement error. A similar phenomenon applies 
for Sobol' indices. In particular, no square GSI or sum of squares will yield 
for |u| < d, as the next proposition shows. 

Proposition 1. The coefficient of erf, in Ylr=i //Eu ^r,uf{x u :z- u )) 2 dx dz 
Proof. It is enough to show this for R = 1 with Xi u = A„. First, 



(J2 \uf(x u :z_ u )) 2 dx dz = A T 9A. 



(12) 



Next, the only elements of containing a 2 -, are the diagonal ones, equal to a 1 . 
Therefore the coefficient of erf, in (fT2l is ^ u \\. □ 

For a square or sum of squares to be free of cr|, it is necessary to have 
X^=i ^ru — 0- That in turn requires all the A rjM to vanish, leading to the 
degenerate case Q = 0. As a result, we cannot get an unbiased sum of squares 
for any GSI that does not include cnk. In particular, cannot have an unbiased 
sum of squares estimate for |u| < d. 



5 Specific variance components 

For any w C T> the variance component for w is given in ([5|) as an alternating 
sum of 2H lower Sobol' indices. It can thus be estimated by a simple GSI, 

f(x) \vf(xv:*-v) (13) 

vC.w 

where = (-l)^-"!. The cost of this simple GSI is C = 2^ + l\ w \ <d . If 
w = T>, then f(x) appears twice, but otherwise it is only used once. The GSI 
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can also be estimated by some bilinear GSIs using fewer function evaluations as 
we show here. 

We begin by noting that for u, v C w, 



NXOR(u, v + w c ) = (XOR(u, v) + w c f = NXOR(u, v)Hw 

and so 

Qu,v+w c = H + ZNXOR(«,u)ntu 

does not involve any of the variables Xj for j £ w. Note especially that 
NXOR(u, v) itself contains all of w c and is not helpful in estimating cr^ when 
\w\ < d. 

To illustrate, suppose that w — {1,2,3}. Let u,v C w. Then we can work 
out a 2 x 8 submatrix of the Sobol' matrix using 



NXOR(«,ti+to c ) 1 2 3 12 13 23 123 

I" 123 23 13 12 3 2 1 

1 23 123 3 2 13 12 1 



(14) 



where we omit braces and commas from the set notation. 

It follows now that we can use a non-simple bilinear GSI, A T G>7, where A 
and 7 are given by 






1 


2 


3 


12 


13 


23 


123 


1 


-1 




















1 





-1 


-1 








1 






(15) 



with u and v — w c given along the top labels in (|15p while the remaining 2 d — 8 
elements of A and 7 are all zero. Specifically, the expected value of 

]T (-l) H+H f(Xu-Z-u)f(x v+w *:z v o- w ) (16) 

uC{l} vC{2,3} 

is cr^ 2 3 i . While equation (p~3|) for w = {1, 2, 3} requires 9 function evaluations 
per (x, z) pair, equation (|16|) only requires 6 function evaluations. For \w\ < d, 
the u = and v — evaluations are different due to the presence of w c , so no 
evaluations are common to both the A and 7 expressions. There are also two 
variants of (fT5"|) that single out variables 2 and 3 respectively, analogously to 
the way that (fTBl) treats variable 1. 

In general, bilinear GSIs let us estimate using 2 k + 2H~ fc function eval- 
uations per (x, z) pair for integer 1 ^ k < \w\ instead of the 2> w > evaluations 
that a simple GSI requires. 

Theorem 3. Let w be a nonempty subset of D for d ^ 1. Let f € L 2 [0, l] d . 
Choose w\ C w and put w-i — w — w\. Then 

°™= E E jj f( x m-z- ui )f(x u2+w .:z u ^ w )dxdz. (17) 
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Proof. Because the right hand side of (ITT)) is a contrast, we can assume that 
H = 0. ThenNXOR(ui,M 2 +w c ) = NXOR(ui, u 2 )nw = w-{u 1 +u 2 ). Therefore 

E E (-i) |uit+M e^ 2+ ^ - E E (-i) M+M ^- ( u 1+tl2) 

= E E (-i) hl " llll+|lU2 ^ 2l ^ 1+U2 

after a change of variable from Uj to lUj — Uj for j = 1,2. We may write the 
above as 

£ £ (.^K-uxl+l^-x*! a l (is) 

Consider the set u C T>. The coefficient of er 2 , in (|T5| is if v n w c 7^ 0. 
Otherwise, we may write v = v\ + t>2 where Vj C tOj, j = 1,2. Then the 
coefficient of o\ 2 in (fTSl) is 



E E (-i) K ~ uil+l " 2| i«c„ 1+U2 

_ (_i)|wi-«i| (_i)l«J2-«2i_ 

Ui :ui Cui Cu;i -^2 :t>2 ^"^2 £^2 

These alternating sums over Uj with Uj C Uj C u^- equal 1 if = ujj but 
otherwise they are zero. Therefore the coefficient of cr 2 , in (|18j) is 1 if i> = w and 
is otherwise. □ 

We can use Theorem [3] to get a bilinear (but not square) estimator of er|, = 
Y 2 -,. A similar argument to that in Theorem [3] yields a bilinear estimator of 
superset importance Tj, for a general set u;. 

Theorem 4. Let w be a nonempty subset of T> for d ^ 1. Let / <= L 2 [0, l] d . 
Choose w\ C w and put u>2 = w — wi . T/ien 

T ^= E E (-i) M+M e^ +tll!l „ C+U2 . (19) 

Proof. Because ai/0, the estimate is a contrast and so we may suppose /i = 0. 
For disjoint u\,U2 C to, NXOR(w c + ui,u> c + M2) = 2? — wi — 1/2, and so the 
right side of (fl9|) equals 

E E (-i) M+M ^- M1 - U2 -EE (-i) |Mil+l " 21 E ^ 

ttiCiui U2QW2 U1C-W1 U2C.W2 vC.T) — ui—U2 

= E-' E E (-i) M+l " 2 W ul _ U2 . 

U til QlOi ti2 Cl02 

Now write v — (v n u; c ) + v 1 + i>2 with C toj, j = 1,2. The coefficient of <7% is 

( — l)l Ul l + l U2 ' l Ul nv 1 ^0^u 2 nv2^0 ■ 

U\ C-W\ U2 ClL>2 
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Now 



(~ l)' ui 'luini)i=0 — 



11 Cwi 



E 



(_!)!«! 



which vanishes unless w% — v% and otherwise equals 1. Therefore the coefficient 
of a 2 is 1 if v D w and is otherwise. □ 

The cost of the estimator ([Jj| is C = 2\ Wl \ +2l™ 2 l - 1, because the evaluation 
f(x w c-.z w ) can be used for both u\ = and U2 = 0. 



6 GSIs with 0(d) function evaluations per pair 

Some problems, like computing mean dimension, can be solved with 0(d) differ- 
ent integrals instead of the 0(2 d ) required to estimate all ANOVA components. 
In this section we enumerate what can be estimated by certain GSIs based on 
only 0(d) carefully chosen function evaluations per (xi,Zi) pair. 



6.1 Cardinality restricted GSIs 

One way to reduce function evaluations to 0(d) is to consider only subsets u 
and v with cardinality 0, 1, d — 1, or d. We suppose for d ^ 2 that j and k are 
distinct elements of V. Letting j and k substitute for {j} and {k} respectively 
we can enumerate NXOR(u, v) for all of these subsets as follows: 



NXOR 


j 


k 


-j 


-k 


V 





' v -j 


-k 


3 


k 





j 


—j v 







{j,k} 


3 


-j 


j 


{j,k} 


V 


~{.hk} 


-j 


V 


j 


k 


-j 


-k 


V 



and hence the accessible elements of the Sobol' matrix are: 








3 





' a 2 




j 




a 2 


-j 







V 





4 



ft 



cr 



-UM 

z -{jM 



V 





(20) 



Using (|2"0)l we can construct estimates of ^ u \u\a 2 — J2j T \j} an d S| u |=i cr 



J2j Lux at cost C = d- 



4< 



1. Simple GSI estimates are available using u — and 



cither v — {j} or v = —{j} for j = 1, 



—2 



, d. More interestingly, it is possible to 
2 n and T 2 . } for j = 1, . . . , d, 
I as was first shown bv lSaltellil (|2002i Theorem 1) 



y fe j along wit h all o 



compute all 1)/2 indices r 

at total cost C = d + 2 as was first shown by ISaltellil (|2002l Theorem 1). Given 
C = 2d +2 evaluations one can also compute all of the t 2 j k ^ indices by pairing 



up u = {j} and v = -{k} (jSaltelli l2002l Theorem 2 
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For the remainder of this section we present some contrast estimators. The 
estimate 



1 

2^Y,(f( x ) ~ f( x {-i} :z {j})) ■ 



is both a contrast and a sum of squares. It has expected value ^2 U \u\a^ and 
cost C = d + 1. 

Next, to estimate X^uIimM^u by a contrast using d + 2 function evaluations 
per (xi, Zi) pair, let 

'i, H = i 

— d, Itil = d. 



A,;, 



Then the contrast X u f(x u :z- U )f(z) has expected value 

E4 = E c7 m = E a «- 

i=i j'=i |u|=i 

The total of second order interactions can be estimated with a contrast at 
cost C = 2d + 2. Taking 



i, M = i 

A u = < -d, |u| = and j v = < 
0, else, 



1, \v\ = d-l 

-(d-2), \v\=d 
0, else 



we get a contrast with 

d d d d 

J=l fe=l fe=l i=l 
= 2j2zl-(2d-2)J2zl = 2j2-l 

|u|=2 |u|=l |u|=2 

Thus A T §7/2 estimates X)|«|=2 at cost C = 2d + 2. 
Next, taking 



A u = < 



1, M = l 



= 1 



0, else, 
we get a contrast with 



:0 and i v ={-(d-2), \v\=0 

else 



d d 



a t ©7 = EEi'm 1 ^^ + E ° 2 d E4> - ( d 2 ) E4> + < rf - 2 ) ff2 

j=l fe=l i=l fe=l J=l 



^ - E X>W^* - 2d(rf - Ik 2 + 2(d - 1)^4 + d(d - 2) ( r 2 

j=l fe=l i=l 
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j=l j — 1 k—1 u 

using Theorem[TJ Therefore E(A T §7/2) = J2 U H 2<J l, at cost C = d+l. 
6.2 Consecutive index GSIs 

A second way to reduce function evaluations to O(d) is to consider only subsets 
u and v of the form {1,2, ... ,j} and {j + 1, . . . ,d}. We write these as (0, j] 
and (j, d] respectively. If f(x) is the result of a process evolving in discrete 
time then (0, j] represents the effects of inputs up to time j and (j, d] represents 
those after time j. A small value of t 2 then means that the first j inputs are 
nearly forgotten while a large value for r? ^ means the initial conditions have 
a lasting effect. 

We suppose for d ^ 2 that 1 ^ j < k ^ d. Once again, we can enumerate 
NXOR(u, v) for all of the subsets of interest: 



NXOR 


S3 


(0,j] 


(0,fc] 




{k,d] 


V 





V 




(k,d] 


(o,i] 


(0,k] 





(0,3] 




V 









(o,j 


U-d] 


(o,j] 





ti,k] 






CM 


V 





(o,i] 


(0,fc] 




(k,d] 


2? 



(21) 



and hence the accessible elements of the Sobol' matrix are n 2 + for the sets 
u in the array above. 

The same strategies used on singletons and their complements can be applied 
to consecutive indices. They yield interesting quantities related to mean dimen- 
sion in the truncation sense. To describe them, we write \u\ = min{j | j E u} 
and \u~\ = max{j | j € u} for the least and greatest indices in the non-empty 
set u. 

Proposition 2. Let / G L 2 [0, l] d have variance components a 2 . Then 

d-1 

E( (M,f - ©0,x>) = ^2{d- \u\)o 2 u , and, 

j = l uCT> 
d-1 

E( e (^i^ - e w) = E (N - 

j=l uCT> 

Proof. Since these are contrasts, we may suppose that fi — 0. Then, using (|2ip 

E e ("^ = E^] = ( d -D- 2 -E^]- 

3=1 J'=l 3=1 
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Next, 6 0i x> = 0, and 

d-l d-l 

^2 T U,d] = CT « 1 {i+i,-,<i}nu#0 = / J(N ^ i) 0- ^- 

3=1 tiCP j=l uCX> 

Combining these yields the first result. The second is similar. □ 

Using Proposition [2] we can obtain an estimate of [it] <7 2 /<7 2 , the mean 
dimension of / in the truncation sense. We also obtain a contrast 

d-l 

^(e^~e (Oj] , c -e (W + 2e ,p) = ^(M - [«JK 

3=1 u 

which measures the extent to which indices at distant time lags contribute im- 
portant interactions. 

We can also construct GSIs based on pairs of segments. For example, 

d 1 d 1) d—1 d 

zl l2 e (o,j],(M 2 ^ 2 = zJ cr «X/ 1 «c(j\/c] 

j=o k=j+i u j=o fc=j+i 

= E°2M(d-M+i)- 



7 Bias corrected GSIs 

When we are interested in estimating a linear combination of variance compo- 
nents, then the corresponding GSI is a contrast. Sometimes estimating a con- 
trast requires an additional function evaluation per (cc^Zj) pair. For instance 
the unbiased estimator ^ of t 2 requires three functio n evaluations per p air 
compared to the two required by the biased estimator of Ijanon et al.l ( 20121) . 

Proposition [3] supplies a bias-corrected version of Janon et al.'s (2011) esti- 
mator of t 2 using only two function evaluations per (xi, Zi) pair. 

Proposition 3. Let f G L 2 [0,l] d and suppose that Xi,Zi ~ U(0, l) d for i = 
1, . . . , n where n ^ 2. Let y i = x i u \Zi _ u and define 



i " i 

-£/(*<), £' = -£/(%), 

i=l i=l 

£(/( Xi )-A) 2 , and S ' 2 = — - A') 2 
i=l i=l 



Then E(r 2 ) = r 2 where 

2n /lA„ , „, , //t + £'\ 2 s 2 + s' 2 
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Proof. First E(/(xi)/( V< )) = M 2 + !«• Next 

E((/t + /t') 2 ) = 4/i 2 + Var(/i) + Var(A') + 2Cov(A, A') 

= 4/ i 2 + 2— +2^. 

n 7i 

Finally E(s 2 ) = E(s' 2 ) = a 2 . Putting these together yields the result. □ 

More generally, suppose that E(tr(Sl T 0)) contains a contribution of /i 2 l T J71 
which is nonzero if fl is not a contrast. Then a bias correction is available for 

tr(n T e). 

Proposition 4. Let f e L 2 [0, l] d and suppose that Xi,Zi ~ U(0, l) d for i = 
1, . . . , n for n 2. For u C P define 

^ 71 1 71 

Au = - J! f( X hU :Z h-u)> and S u = _ 1 / l(/(jgj,u'^i,-u) - An) 2 - 
i=l i=l 

2n v^V^r, /A /A« + Ai>\ 2 , s « + s « 2. 



2n 



is an unbiased estimate of J2 U ^ra(9u» ~ A* 2 )- 

Proof. This follows by applying Proposition [3] term by term. □ 

The computational burden for the unbiased estimator in Proposition|4]is not 
much greater than that for the possibly biased estimator tr(f2 T 0). It requires 
no additional function evaluations. The quantities fi u and s 2 need only be 
computed for sets u C T> for which VL UV or Q vu is nonzero for some v. If is a 
sum of bilinear estimators then the — J2 U ^u«A«Au/2 cross terms also have 
that property. 

The bias correction in estimator (|22|) complicates calculation of confidence 
intervals for tr(i7 T 9). Jackknife or bootstrap methods will work but confidence 
intervals for contrasts are much simpler because the estimators are simple aver- 
ages. 



8 Comparisons 



There is a 2 2d -dimensional space of GSIs but only a 2 d — 1-dimensional space 
of linear combinations of variance components to estimate. As a result there is 
more than one way to estimate a desired linear combination of variance compo- 
nents. 

As a case in point the Sobol' i ndex r 2 can be estim ated by either the original 



method or by the contrast Q. iJanon et al.l (|2012f ) prove that their estimate 



of A improves on the simpler one and establish asymptotic efficiency for their 
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estimator within a class of methods based on exchangeability, but that class 
does not include the contrast. Similarly, inspecting the Sobol' matrix yields 
at least four ways to estimate the variance component a% 23 j, and superset 
importance can be estimated via a square or a bilinear term. 

Here we consider some theoretical aspects of the comparison, but they do 
not lead to unambiguous choices. Next we consider a small set of empirical 
investigations. 



8.1 Minimum variance estimation 

Ideally we would like to choose f2 to minimize the variance of the sample GSI. 
But, the variance of a GSI depends on fourth moments of AN OVA contribu- 
tions which are ordinarily unknown and harder to estimate than the variance 
components themselves. 

The same issue comes up in the estimation of variance components, where 
MINQE (minimum norm quadratic estimation) estimators were proposed in a 
seri es of papers by C. R . Rao in the 1970s. For a comprehensive treatment 
see iRao and Kleffej (|l988h who present MINQUE and MINQIE versions using 



unbiasedness or invariance as constraints. The idea in MINQUE estimation is 
to minimize a convenient quadratic norm as a proxy for the variance of the 
estimator. 

The GSI context involves variance components for crossed random effects 
models with interactions of all orders. Even the two way crossed random effects 
model with an interaction is complicate d enough tha t no closed form estimator 



appears to be known for that case. See iKleffe 



gh tha t 
(1980). 



We can however generalize the idea behind MINQE estimators to the GSI 
setting. Writing 

Var(tr(ft T 9)) = £ X! X! X! ^JWCov^, §„,„,) 

nCPtiCDu'CDo'CP 

we can obtain the upper bound 



Var(tr(f> T §)) ^ f X X M™| 2 ) max Var(9 ! „) , 

leading to a proxy measure 

V(tl) = X l^™! 2 = tr(tt T ft). (23) 



uCVvCV 



Using the proxy for variance suggests choosing the estimator which min- 
imizes C(fi) x F(O). The contrast estimator (gj) of zt has C(fi) x V{Q) = 
3x2 = 6 while the original Sobol' estimator has C(fl) x V(O) =2x1 = 2. 
The estimators ([TBf and p^]) for 2 3 , both have V(Q) — 8. The former has 
cost C(f2) = 9, while the latter costs C(fi) =6. As a result, the proxy argu- 
ments support the original Sobol' estimator and the alternative estimator (fT6t 

fOT CT {1,2,3}' 
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8.2 Test cases 

To compare some estimators we use test functions of product form: 

d 

f(*)=tl(to+r j gi(x j )) (24) 
i=i 

where each gj satisfies 

g{x) dx = 0, / g(x) 2 dx = 1, and / g(x) 4 dx < oo. 
o Jo Jo 

The third condition ensures that all GSIs have finite variance, while the first 
two allow us to write the variance components of / as 



2 = jn je „r?xn^ u /4 h>o 

else, 



(7,, = 



along with /i = rX/_i Mj- 

We will compare Monte Carlo estimates and so smoothness or otherwise of 
<7j(-) play no role. Only fij, Tj and the third and fourth moments of g play 
a role. Monte Carlo estimation is suitable when / is inexpensive to evaluate, 
like surrogate functions in computer experiments. For our examples we take 
g 3 (x) = VT2(x - 1/2) for all j. 

For an example function of non-product form, we take the minimum, 

f(x) = min Xj. 



Liu and Owenl ( 2006h show that 



(d+l) 2 (2d-| U | + 2)' 
for this function. Taking u — V, gives a 2 = d(d + l)~ 2 {d + 2) _1 . 



8.3 Estimation of cr^ 2 3 i 

We considered both simple and bilinear estimators of cx^ 2 3 j in Section [S] The 
simple estimator requires 9 function evaluations per (ce, z) pair, while three 
different bilinear ones each require only 6. 

For a function of product form, all four of these estimators yield the same 
answer for any specific set of (xi,Zi) pairs. As a result the bilinear formulas 
dominate the simple one for product functions. 

For the minimum function, with d — 5 we find that by symmetry, 

°l = lfi l2)3} - 3^1,2} + 3lf 1} = 5^5 = 1-68 x 1CT 4 . 
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Estimator 


Simple 


Bilin.{l} 


Bilin.{2} 


Bilin.{3} 


Mean 


1.74 x 10" 4 


1.72 x 10" 4 


1.68 x 10" 4 


1.70 x 10" 4 


Standard error 


1.05 x 10~ 5 


5.69 x 10~ 6 


5.71 x 10~ 6 


5.67 x 10~ 6 



Table 1: Estimated mean and corresponding standard error for three estimators 
of cr?! 2 3} f° r f( x ) = niinisgj^ Xj when x ~ U(0, l) 5 . The Bilin.jl} estimator 
is from equation <jTH]) , and the other Bilinear estimators are defined analogously. 



Because we are interested in comparing the variance of estimators of a vari- 
ance, a larger sample is warranted than if we were simply estimating a variance 
component. Based on 1,000,000 function evaluations we find the estimated 
means and standard errors are given in Table [1] We see that the bilinear esti- 
mators give about half the standard error of the simple estimator, corresponding 
to about (1.05/.571) 2 x 9/6 = 5.1 times the statistical efficiency. 

8.4 Estimation of 2 | 

We consider two estimators of r 2 . The estimator (j4|) is a bilinear contrast, 
aver aging f(x)(.f(x u :z u ) — f(z)). The estimator ([3J using the estimator of /t 



from Janon et al.l (2012J) is a modification of Sobol's original simple estimator 



based on averaging f(x)f(x u :Z- u ). The bias correction of Section [7] makes an 
asymptotically negligible difference, so we do not consider it here. 

The contrast estimator requires 3 function evaluations per {x, z) pair, while 
Sobol's only requires 2. Both estimators make an adjustment to compensate for 
the bias /j?. Estimator ([3]) subtracts an estimate jx 2 based on combining all 2n 
function evaluations, the square of the most natural way to estimate \i from the 
available data. Estimator ((4]) subtracts (1/n 2 ) J2i Yle f( x i)f( x i,u'- z i,-u), which 
may be advantageous when the difference f (x u :Z- u )— f(z) involves considerable 
cancellation, as it might if x u is unimportant. Thus we might expect (|4j) to be 
better when r 2 is small. We compare the estimators on a product function, 
looking at r 2 for three subsets u of size 2 and varying importance. 

For the product function with d = 6, r = (1, 1, 1/2, 1/2, 1/4, 1/4), all fij = 1 
for j = 1, . . . , 6, and gj(xj) = \/Y2{xj — 1/2), we may compute j3 1 2 } = 3 = 
0.50cr 2 , t 2 3 4} = 1.56 = 0.093ct 2 , and r 2 5 6} = 0.13 = 0.021ct 2 . 

Results from R = 10,000 trials with n = 10,000 (a^, Zi) pairs each, are shown 
in Table [2] The efficiency of the contrast estimator compared to the simple one 
ranges from about 0.5 to about 2.5 in this example, depending on the size of 
the eff ect being estimated , with the contrast being better for the small quantity 



Tr 5 6 i . ISobol et al.l (|2007l ) also report superiority of the contrast estimator on a 



Neither estimator is always more efficient than the other, hence no proxy 
based solely on f2 can reliably predict which of these is better for a specific 
problem. 
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n = 10,000 


{1,2} 


{3,4} 


{5,6} 




Cont. 


Simp. 


Cont. 


Simp. 


Cont. 


Simp. 


True 
Avg. 
Bias 


3.0000 
3.0002 
0.0002 


3.0000 
3.0002 
0.0002 


0.5625 
0.5624 
-0.0005 


0.5625 
0.5628 
0.0003 


0.1289 
0.1291 
0.0002 


0.1289 
0.1294 
0.0005 


S.Dev 

Neg 


0.1325 


0.1186 


0.0800 


0.0998 


0.0378 
0.0001 


0.0737 
0.0336 


Eff. 


0.53 


1.04 


2.54 



Table 2: This table compares a simple estimator versus a contrast for t 2 with 
n = 10,000. The sets compared are u = {1,2}, {3,4}, and {5,6} and / is 
the product function described in the text. The rows give the true values of 
r 2 , and for 10,000 replicates, the (rounded) sample values of their average, bias, 
standard deviation and proportion negative. The last line is estimated efficiency 
of the contrast, (2/3) times the ratio of the standard deviations squared. 

The bias correction from Section [7] makes little difference here because for 
n = 10,000 there is very little bias to correct. It does make a difference when 
n = 100 (data not shown) but at such small sample sizes the standard deviation 
of can be comparable to or larger than r 2 itself for this function. 

8.5 Estimation of 2 3 4 j 

Here we compare two estimates of T 2 ^ 234}, the square (|11[) and the bilinear es- 
timator dT9l from TheoremSJ For a product function, = Ylj ew rj Iljgu; 
tJ). Squares have an advantage estimating small GSIs so we consider one small 
and one large (for a four way interaction) T 2 . 

For d = 8, t = c(4,4,3,3,2,2, l,l)/4 and all fij = 1 we find that T 2 {1 2 3 4} = 
0.558 = 0.0334ct 2 and T 2 5 6 7 8} = 0.00238 = 0.000147cr 2 . The bilinear esti- 
mate (fT9"|) based on w\ — {1,2} and W2 — {3,4} for T{ li2 ,3,4} (respectively 
Wi = {5, 6} and W2 = {7, 8} for T{ 5 6 .7,8}) requires C — 7 function evaluations, 
while the square (jTTJ) requires C = 16. From Table[3]we see that the square has 
an advantage that more than compensates for using a larger number of function 
evaluations and the advantage is overwhelming for the smaller effect. 

The outlook for the bilinear estimator of is pessimistic. Its cost advantage 
grows with \w\; for \w\ = 20 it has cost 1023 compared to 2 20 for the square. 
But T 2 ^ for such a large w will often be so small that the variance advantage 
from using a square will be extreme. 
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{1,2,3,4} {5,6,7,8} 



Bilinear 
Square 



35.07 
6.04 



4.019 
0.051 



Efficiency 



14.7 



2,710 



Table 3: Standard errors for estimation of Tj, by the bilinear estimate and 
a square as described in the text. The estimated standard errors based on 
n = 1,000,000 replicates are 10~ 3 times the values shown. The relative efficiency 
of the square is 7/16 times the squared ratio of standard deviations. 

9 Conclusions 

We have generalized Sobol' indices to estimators of arbitrary linear combina- 
tions of variance components. Sometimes there are multiple ways to estimate 
a generalized Sobol' index with important efficiency differences. Square GSIs 
where available are very effective. When no square or sum of squares is available 
a bilinear or low rank GSI can at least save some function evaluations. Contrasts 
are simpler to work than other GSIs, because they avoid bias corrections. 
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